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ANALYSIS OF THE TWO-REGIME METHOD ON SQUARE MESHES 

MARK B. FLEGG*, S. JONATHAN CHAPMAN* , LIKUN ZHENGt, AND 

RADEK ERBAN* 

Abstract. The two-regime method (TRM) has been recently developed for optimizing stochas- 
tic reaction-diffusion simulations |11| . It is a multiscale (hybrid) algorithm which uses stochastic 
reaction-diffusion models with different levels of detail in different parts of the computational do- 
main. The coupling condition on the interface between different modelling regimes of the TRM was 

£/") previously derived for one-dimensional models. In this paper, the TRM is generalized to higher di- 

r" H mensional reaction-diffusion systems. Coupling Brownian dynamics models with compartment-based 

models on regular (square) two-dimensional lattices is studied in detail. In this case, the interface 

£vj between different modelling regimes contain either flat parts or right-angled corners. Both cases are 

studied in the paper. For flat interfaces, it is shown that the one-dimensional theory can be used 
along the line perpendicular to the TRM interface. In the direction tangential to the interface, two 
choices of the TRM parameters are presented. Their applicability depends on the compartment size 

^J>( and the time step used in the molecular-based regime. The two-dimensional generalization of the 

1^^ TRM is also discussed in the case of corners. 
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^*^ 1. Introduction. There are two common approaches to stochastic reaction- 

diffusion modelling: (i) compartment-based models; and (ii) molecular-based models 
[HIE]- Molecular- based models provide a higher level of detail, but they are often 
more computationally intensive than compartment-based models. In some applica- 
tions, microscopic detail is only required in a relatively small region, for example, 
close to the cellular membrane or a particular organelle [T^J IS] • Such problems are 
best simulated by a hybrid method which uses a detailed modelling approach in lo- 
calized regions of particular interest (in which accuracy and microscopic detail is im- 
portant) and a less detailed model in other regions in which accuracy may be traded 
for simulation efficiency. To apply this general idea to stochastic reaction-diffusion 

l/~j modelling, one has to introduce a suitable boundary condition between different mod- 

I* elling regimes. In [llj . we derived the appropriate boundary condition for coupling 

one-dimensional compartment-based and molecular-based models. We developed the 
two-regime method (TRM) which has the accuracy of the detailed molecular-based 
approach (in the region where it is required), but benefits from the efficiency of a 
^. less detailed (coarser) compartment-based model in other parts of the computational 

domain. In this paper, we generalize the TRM to higher dimensional simulations. 

In the remainder of this first section, we introduce the notation which is used 
during the rest of this manuscript. We summarize both compartment-based (Section 
|l.ip and molecular-based modelling (Section |1.2| . Then we introduce the TRM in 
Section |1.3| Our main results are presented and derived in Section [2] In Section 
|3j we demonstrate the applicability of the presented theory using several illustrative 
numerical examples. Finally, we discuss other multiscale (hybrid) stochastic reaction- 
diffusion approaches from the literature. These methods vary in implementation and 
applicability and we put them into context against the TRM in Section HI 
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1.1. Compartment-based modelling. Mesoscale compartment-based mod- 
elling of reaction-diffusion processes begins by partitioning the domain Qc m t° com- 
partments (open sets) Cj, j = 1, . . . , K, such that the compartments do not overlap 
and they cover the whole domain ilc (i.e. uf =1 Cj = ilc and Ci n Cj = 0, for i ^ j, 
where overbars denote the closure of the corresponding set) . Assuming that there are 
M chemical species Zi, i = 1, . . . , M, the state of the system is completely defined 
by the copy numbers Mi j £ No of molecules for chemical species Zi found in the 
compartment Cj, i = 1, . . . , M, j = 1, . . . , K . In what follows, symbol No denotes 
the set of nonnegative integers, i.e. No = {0, 1,2,3,... }. The simulation of reaction 
and diffusion of the molecules in the system is usually implemented by event-driven 
algorithms, which include the Gillespie algorithm [17], the Next Reaction Method 
[TB] or the Next Subvolume Method [20 . They have been implemented in several 
open-source software packages including MesoRD [2D], URDME [5], STEPS [3JJ and 
SmartCell [2J. In this paper, we will use a derivative of the Next Reaction Method 
from Gibson and Bruck 16 . 

Event-driven algorithms require the calculation of event propensities for a par- 
ticular system state [17]. An event propensity ag is the rate (per unit time) for an 
event £ to occur that changes the state of the system. Events in compartment-based 
reaction-diffusion processes may include: reaction events (in which chemical molecules 
of some species change into molecules of other species, or are just introduced or re- 
moved from the system), diffusive events (in which molecules of a chemical species 
jump from one compartment to an adjacent compartment) or boundary events (in 
which molecules are absorbed by, react with or reflect from a domain boundary). A 
putative time tg for each event £ can be found given some current time t using 

t e =t+—hx(—Y (i.i) 

where T£ are uniformly distributed random numbers between and 1 chosen separately 
for each occurence of each event. The next event which takes place in the system is 
determined by finding which event corresponds to time mim; t$ where the minimum 
is taken over the set of all possible events [IB]. The state is changed to reflect the 
occurence of the event and the current time is then updated to t := miiigtg. The 
current event might also change propensities of some related events. The putative 
times for these events must therefore be scaled to reflect the change in propensity. 
That is, 



tT W --=t+-L,(4 d -t), (1.2) 
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where a^ ld and a^ cw (resp. t| ld and ££ cw ) are the propensities (resp. putative times) 
for event £ before and after the current event takes place [IB]. A putative time 



for the next occurence of the current event must be resampled using formula (1.1). 
The simulation is then constructed by a series of successive events over time, in each 
instance, defined by the most imminent event and performing the state change that 
defines that event. 

1.1.1. Reaction events. One of the main assumptions of compartment mod- 
elling of stochastic reaction-diffusion models is that each compartment is small enough 
that it may be considered well mixed [S] . Reactions are modelled in each compartment 
by defining the propensity for reaction in each compartment. Consider the reaction 



Table 1.1 
Examples of propensities and effects of reactions in compartments. To simplify this table, we 
omit species with zero coefficients in reaction and product complexes which were included in the 
general form jl.3| . 



example reaction TZ 


tt-Rj 


changes of the state vector 


Zi + Zk — > Zi 




A/ij changes to A/ij — 1 
Nk.j changes to A4,j — 1 
A/jj changes to A/jj + 1 


2Z % A Zi 


kM.j(Mj - 1) 


A/ij changes to A/ij — 2 
A/"; j changes to A/jj + 1 


z % A z k + z t 


nAfi,j 


Ni,j changes to A/ij — 1 
A4,j changes to A4,j + 1 
Ml j changes to A/jj + 1 


0Az, 


kVj 


A/jj changes to A/jj + 1 



TZ in compartment Cj given by the general form 
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where /% 6 No (resp. 7* G No) are the numbers of molecules of chemical Zi, i = 
1,2,..., M, that are required as reactants (resp. products) of the chemical reaction 
TZ and k is the reaction rate. We define the notation for this event £ = (7Z,j). In 
realizing this reaction event in the j-th compartment, the number of molecules A/jj, 
i = l,2,..., M, j = 1, 2, . . . , K , change by the corresponding stoichiometric coefficient 
vi = 7,- — /3j. Considering mass action chemical kinetics, the propensity for this event 
to occur depends on the number of available reactants in the compartment Cj. In 3D, 
we can postulate this dependence in the following form [HI [7] 



rl-EiA 
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(1.4) 



where Vj is the volume of the compartment Cj . Table |1.1| shows the examples of 
propensity an.j for some simple reactions TZ in compartment Cj and the subsequent 
changes to the state of the system that occur as a result of the reaction event. 



The main focus of this paper is on 2D simulations. In this case, the formula ( 1.4 1 



has to be slightly modified by replacing the compartment volume Vj by its area. Since 
the propensity u-r.j is dimensionless, some rate constants have different physical units 
in 2D than in 3D. In applications, 2D simulations in the domain ile C M 2 are often 
viewed as a model of a real 3D domain flc X (0, w) where the domain width w is so 
small that the spatial distribution along the third axis can be neglected. Then the 



formula (1.4) can be applied with the standard interpretation (i.e. physical units) of 
the reaction rates. The area of each 2D compartment is multiplied by w to get the 



corresponding volume Vj in (1.4). 

1.1.2. Diffusion events. Diffusion in compartment-based models of reaction- 
diffusion processes is defined by the stochastic jumping of molecules of chemical species 



Zi, i — 1,2,..., M, between any two adjacent compartments from Cj to Ck, j, k = 
1,2, ... ,K. We shall define the notation for this diffusive event to be £ = {T>, i, j, k). 
The propensity for a diffusive event is given by 

<*D,ij,k = 1j,kDiNij, (1.5) 

where qj t k is dependent on the morphology and relative positions of the compartments 
Cj and Ck, and Di is the diffusion constant for chemical Zi. For adjacent square or 
cubic compartments Cj and Ck of length h on a regular lattice, we have qj t k = 1/h 2 , 
i.e. 

ot-D,i,j,k = T^M,jj (1-6) 

and qj^k = if Cj and Ck do not share a common side. For more irregular com- 
partment shapes, qj t k can be determined by the finite element discretisation of the 
diffusion equation on a lattice whose vertices are at the centers of the compartments 
[5]. During the diffusive event {V, i,j, k), the state of the system is changed to reflect 
the movement of one molecule of Zi from Cj to Ck, i.e. Afij changes to A/ij — 1 and 
Mi,k changes to A/i,fc + 1- 

1.1.3. Boundary events. Boundary events are caused by the diffusion of mole- 
cules into the boundaries of the domain dilc- It is no surprise then that boundary 
events are linked closely to diffusion events. Consider the compartment Cj which is 
adjacent to dQc ■ A diffusive event of a molecule of chemical Zi which would ordinarily 
result in a jump from Cj to a compartment on the other side of dVLc actually results in 
an interaction of the molecule with the boundary. Such an interaction usually results 
in one of two outcomes either: the molecule is absorbed; or the molecule is reflected 
[S] . Molecules attached to the surface can be released back into the solution [3J [23] . 
However, in this paper, we will consider reflective boundary conditions on all external 
boundaries for simplicity. 

1.2. Molecular-based modelling. Molecular-based approaches to reaction- 
diffusion modelling are characterized by the prescription of exact coordinates in space 
for each molecule of each chemical species Z i: i — 1,2, ... ,M, in the continuous 
domain Qm C R w , where N — 1,2,3. The trajectory of large molecules (such as 
proteins) are computed using Brownian dynamics 4, 14, 25J. A random displacement 
of each molecule every timestep At models the effect of solvent molecules on the large 
molecules of interest without the need to simulate each solvent molecule individually. 
We will denote the j'-th molecule of chemical species Zi as Zf . Given a time step At 
the position Xjj(t + At) of molecule Zf at time t + At is computed from its position 
Xjj(t) at time t in the A^-dimcnsional continous space Qm by 

x id (t + At) = x SJ (t) + V2AAtC, (1.7) 

where £ € M. N is a vector of N uncorrelated normally distributed random numbers 
with zero mean and unit variance chosen separately for each molecule. The equation 



(1.7) is a discretization of the standard Brownian motion 

dXij = v / 2AdWy. (1.8) 

Molecular-based models have been implemented in several software packages, includ- 
ing Smoldyn g|, MCell [2S] and GFRD [23] ■ They have been used for modelling 
several biological systems, including the signal transduction in bacterium E. coli [24] 
and the MAPK pathway [27]. 



1.2.1. Chemical reactions. It is relatively straigtforward to implement zero- 
order and first-order chemical reactions in molecular-based models |S]. There is a 
variety of different ways to model bimolecular (second-order) molecular-based reac- 
tions. Consider two molecules Z\ and Z k which can react according to the following 
bimolecular reaction 



Then a suitable probability of reaction per time step is chosen such that the macro- 
scopic rate of reaction between the two chemicals is k [7] . This probability is a function 
of the distance |xjj — Xjy| between reacting molecules. Models postulate that the 
molecules can only react if they are within a specific distance (reaction radius) [4]. 
Care must also be taken when generating the initial positions of products of chem- 
ical reactions. This is especially the case whenever reversible reactions are present 
in the system. If the molecules are not initialized properly then the products may 
unphysically react immediately after being created [23 [5TJ . 

1.2.2. Boundary conditions. Molecules that migrate over domain boundaries 
are treated depending on whether they are reflective, absorbing or reactive boundaries 



[6j. When modelling boundary conditions, one has to take into account that (f.7l 



is only a n ap proximation of the Brownian motion (1.8 1. Let us consider that the 
formula (1.7 1 gives the updated position of the molecule Z\ inside the domain Om, 



i.e. Xjj(i + At) £ £7j\/. Then there is still a nonzero probability that the molecule 



(which follows (1.8)) left the domain Qm during the time (t,t + At) and then returned 
back to Qm. We denote this probability by P m = P m (x ij (i),Xj J (i + At)). Near a 
flat boundary dflM, probability P m takes the analytical form 

p _ ( -distjxj, j(t), dfl M ) dist(x i|J -(t + At),d£l M ) \ n Q , 

where dist(-, OCIm) represents the distance from the boundary 8H,m U- Whilst bound- 
aries are sometimes only implemented at time t + At for molecule Z\ if X4,j(t + At) 



computed by (1.7) is outside Qm, thorough implementation of boundary conditions 
should not only be considered when Xij(t + At) (jz. Qm but also with the probability 
P m if Xi.j(t + At) £ fi m- No closed form solution exists for P m for irregular boundary 
geometries 8Qm- In practice, a curved boundary is usually described locally by a flat 
approximation which increasingly becomes more accurate for small values of At [3] • 

1.3. The two-regime method. Molecular-based techniques are usually pref- 
ered over compartment-based simulation techniques when the concentration of molecu- 
les is low to give a level of microscopic detail that is not achieveable in the mesoscopic 
compartment-based approaches. However, it can be very cumbersome numerically to 
simulate every single molecule and perform probability tests for every pair of molecules 
that have the potential to react if the copy numbers of molecular species are large. In 
such cases, compartment-based approaches are appropriate. They provide a level of 
efficiency that does not require the tracking of each individual molecule. This comes, 
however, at the cost of detail in the simulation. This paper focuses on the TRM which 
provides a way to spatially connect regions that use compartment-based modelling in 
regions that require a mesoscopic level of detail with that of molecular-based mod- 
elling in regions where the concentrations are lower or a high level of detail is required 

El- 



The TRM for simulation of stochastic reaction-diffusion processes is characterized 
by its partition of the computational domain SI C R N , N = 1,2,3, into two non- 
overlapping open subsets Qc and ^m, i.e. fie H fijvf = and fie U fi.M — fi where 
overbars denote the closure of the corresponding set. We denote by I the interface 
between the subdomains fie and fljfi !•©• / = <9fic H <9£1m- Internally, molecules are 
simulated in subdomains fie and fijvf by compartment-based and molecular-based 
approaches which were described in Sections |1.1| and |1.2| respectively. Molecules in 
Qm arc updated at prescribed times separated by At. Meanwhile, molecules in fie 
are updated at the events determined by compartment modelling rules described in 
Section |1.1| The simulation is thus built from a series of time steps that occur at 
each "regular time step" separated by At and each event inside fie- As the time 
updates can be classified by the region that they apply to, updates corresponding to 
compartment-based events are known as C-events (or compartment events) and the 
regular updates separated by At in time are known as M-events. There are several 
possible variants of the TRM [IT] . In this paper, we will analyse the TRM in the form 
which is summarized in Table [L2| In the step (i), we define fie, fi-fcf and the time 
step At. Initial conditions in CIm and fie are implemented in the steps (ii) and (iii), 
respectively. In the step (iv), we also define putative times tM and tc when the next 
M-event and C-event will occur, respectively. Then the TRM repeats steps (v) and 
(vi) until the desired end of the simulation. 

1.3.1. Transition of molecules from fie to VLm . The partition of the domain 
into fie and fijvj is an artificial partitioning that should not interfere with the natural 
diffusion of molecules in the domain. To describe migration of molecules from fie to 
fijvf! we need an expression for the propensity of molecules to jump from compartments 
Cj, adjacent to the interface /, into Om- When a molecule successfully jumps from 
Cj into fijvfj it must be placed with a specific set of coordinates by virtue of the 
modelling approach in VLm- For a regular array of square or cubic compartments in 
fie, we define the propensity of chemical species Zi to jump from Cj into fijvf to 
be some multiple $, j times the natural jumping propensity between neighbouring 
compartments (provided that compartment Cj is adjacent to /, a so-called interfacial 



compartment, otherwise this propensity is equal to 0). That is (compare with (1.6)), 



ocv,i,j,M = $i,j-rjNi,j, (1-10) 

where the index j is a reference to the originating interfacial compartment Cj . Subse- 
quently, the molecule, after being chosen to jump from Cj into Qm is initialized at a 
position x € CIm- We do not restrict this initialization to a specific location but rather 
consider that the initial position is chosen from a probability distribution fij(x). 

1.3.2. Transition of molecules from fiju to fie- In step (vi), a molecule 
originating in VLm is transfered into fie with a probability \P G [0, 1] if the molecule 
interacted with the interface / within the time interval [t,t + At]. It is postulated 
that molecule Z] interacted with / if one of these two conditions is satisfied: 



(a) Xjj(t + At) computed by (1.7) satisfies Xi.j(t + At) <G flc 

(b) Xj ,-(t) G £Im and r < P m where r is uniformly distributed random number 
in (0, 1) and P m = P m (xij(t), x^j(t + At)) was introduced in Section 



1.2.2 



For a straight interface /, the probability P m is given by (1.9). 

If the probability \& is strictly less than 1, then we have to incorporate into the TRM 
that all molecules which satisfy (a) and which are not transported to fie are reflected 



Table 1.2 
The pseudocode of the TRM for stochastic reaction- diffusion simulation. 



(i) Define the subdomains Vtc and VLm and the interface / = dilc H dflm- 
Subdivide O c into compartments Cj, j = 1, . . . ,K. Choose the time step 
At between updates of the molecular-based regime (M-events) in Qm- 

(ii) Specify the initial condition in Qm by placing molecules Zj in U,m at initial 
positions Xjj(O) € fijvf, i = 1,2, ... , M, j = 1,2,... , n(z), where n(i) is the 
initial number of molecules of the i-th chemical species Zi in VLm- 

(iii) Specify the initial condition in Qc by initializing the copy numbers Afij in 
f2c for each chemical species Zi in each compartment Cj, i — 1, 2, . . . , M, 
j = 1, 2, . . . , K. Initialize time as t := 0. 

(iv) Use equation ( |1.1[ ) to calculate ts, the putative times at which all C-events 
£ will take place. Set tu = At and tc = mirif if where the minimum is 
taken over all possible C-events £ . 

(v) If tc < £m, then the next C-event occurs: 

• Update current time t := tc- 

• Change the number of molecules in Q,c to reflect the specific C-event 
that has occured. If this event is one in which a molecule of chemical 
species Zi leaves ftc bound for £Im, then compute its initial position 
in Qm according to the probability distribution /^(x) and remove it 
from the corresponding compartment Cj. 

• Calculate the next putative time for the current C-event by equation, 



(1.1). For all propensity functions ag that are changed as a result of 



the C-event, determine the putative times of the corresponding event 



by equation (1.2) 



• Set t c := mm £ (T £ ). 

(vi) If tM < tc, then the next M-event occurs: 

• Update current time t := tyi- 



Change the locations of all molecules in VLm using equation (1.7). 
Implement boundary conditions at the external boundary dfl^j \ I. 
Initialize all molecules which migrated from flc to f2j\/ since the previ- 
ous M-event at locations computed in the step (v) according to ftj(x). 
Perform all reaction events in Qm- 
Identify all molecules that interact with the interface / from VLm (ex- 



cluding those just initiated) using conditions (a)-(b) from Section 1.3.2 
Move each molecule into the appropriate compartment in flc with prob- 
ability ty. Otherwise, its position is reflected back to J7m- 
For all propensity functions ag that are changed as a result of the 
M-event determine the putative times of the corresponding C-event by 



equation (1.2) 



• Update tM '■= tht + At and, if necessary, set tc '■= mmg(rf ). 
(vii) Repeat steps (v) and (vii) until the desired end of the simulation. 



back to Om- However, this condition is not necessary in ID where it is possible to 
prove that \& = 1 llj. This simplifies the implementation of the TRM in ID. In the 
next section, we summarize the results of the ID theory presented in [llj . Then, in 
Section [2j we analyse the TRM in 2D. 



1.3.3. Summary of analysis of the TRM in ID. In [IT], the TRM is pre- 
sented and analysed in a one dimensional domain fi = (—00, 00) which was divided 
by the interface / = {0} into ilc = (— 00, 0) and Ojvf = (0, 00). The subdomain 
flc was divided into compartments of the same length h, i.e. Cj — {—jh, (1 — j)h), 

j = 1, 2, In this case, there is only one interfacial compartment corresponding to 

j = 1 and coordinates x G Qm are defined as the displacement from the interface I. 
It was obtained that 

•*'A- *"'■ (L11) 



fu(x) = J ^ wfr. , ^ 1 . xefi M , (1.12) 

where erfc(a;) = 2/y/w J exp(— t 2 )dt is the complementary error function. These 
formulae were derived under the assumption DAt ~ ft 2 as ft. — > 0, a condition which 
is used throughout the manuscript. If At gets larger than this then the expected 
jump distance of molecules in the Brownian domain, v 2DAt, becomes greater than 
the compartment size h. Since in the simplest form of the TRM scheme we do not 
allow Brownian particles to jump to interior compartments (particles that cross the 
interface are placed in a boundary compartment), we often assume D At < h 2 (we are 
usually interested in a finer resolution in the Brownian domain than the compartment 
domain). If DAt is chosen much smaller than h 2 then the coupling between domains 
is still accurate in one dimension. However, we will see that choosing \j2DAt <C h 
in higher dimensions can create some numerical artefacts associated with molecules 
diffusing along the interface I. 

2. Main results. In the rest of this paper we wish to derive forms for <j>jj, if? 
and fi t j(x) in domains with dimensions greater than N = 1. Since these parameters, 
the sole requirements for the correct spatial coupling of regions Qc & n d fJjw, are 
independent of the parameters for reactions we do not need to note the chemical 
species in our analysis. We shall therefore drop the index i and relabel Di as D which 
will denote the diffusion constant for the chemical species in question. We therefore 
wish to find the parameters $.,, "J and fj(x) for domains that have more than one 
dimension. We shall limit the analysis to regularly spaced square compartments of 
length h for dimension N = 2, but the results can be easily generalized to regularly 



spaced ./V-dimensional cubic compartments (see Section 2.3). The parameters will be 
derived for a flat interface / and then we will discuss the case where I may have a 
corner. 

2.1. Matching at a flat interface in 2D. The considered geometry is repre- 



sented graphically in Figure 2.1 We present here a derivation to the parameters $j, 
^ and fj (x) for the TRM applied to an infinite two-dimensional domain £1 = R 2 and 
a single flat interface /. To derive these algorithm parameters (and formulate them 
in a reasonably simplified form), we denote by (x,y) the Cartesian coordinates that 
describe domain 57. Without loss of generality we assign Om to the region defined 
by x > and therefore Sic to the region defined by x < (i.e. the interface I is 
the line x = 0) . The compartments Cj are regularly spaced squares of side length h. 
We find it convenient to describe the compartments by two indicies such that each 
compartment Ci.j is assigned to the region described by — (i + l)h < x < —ih and 
jh — h/2 < y < jh + h/2, where i € No and j € Z. In what follows, we will denote 
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Fig. 2.1. Graphical representation of the TRM in 2D for a flat interface I . Compartment-based 
regime Qq is on the left (numbers denote the number of molecules in the corresponding compart- 
ment), molecular- based regime Qm on the right (three illustrative trejectories of individual molecules 
are plotted as red lines). The interface I is plotted as a yellow line. 



by Pij(t)h the probability of a molecule to be in the compartment C ij5 i.e. Pij{t) is 
the (discretized, averaged) probability density in the compartment C^.j. 

Since the analysis of the TRM only depends on the properties of diffusion [TTj . 
we can write the governing equations as the evolution equations for the probability 
density of a single diffusing molecule in 51. The goal of the TRM is to correctly 
approximate its probability density function P(x,y,t) where (x,y) G O and time 
t > 0. Using equations (1.6 1, (1.7), (1.9|, (1.10) and the description of the evolution 



of molecules near the interface / during the TRM, the TRM master equations for both 
the average probability density po. n (t) for a molecule to be in the Co jra compartment 
and the probability density p(x, y, t) in J1m are given by 

A , / (3 + $ )DAt\ rs DAt, ,. ,. , .. 

P0, n {t + At) = I 1 - ^ j^ J P0,n(t) + -p-(P0,n-l(i) +p ,n+l(t) +Pl,n(*)) 



9x7/ /•nh+h/2 



/>nh-\-h/2 /*oo /*oo /*cxd 

/ dy / dy / dx / dx p(x, y, t)K (x + x)K {y - y) , (2.1) 

Jnh-h/2 J-oo JO Jo 



/OO POO 

dy / dx p(x, y, t) [K (x - x) + (1 - 2^)K (x + x)} K (y - y) 
-oo JO 

+ DAtJ2^ j f j (x,y)p , j (t), (2.2) 

where K(x) = {AitDAt^ 1 / 2 exp(— x 2 /(4DAt)) is the distribution of the random dis- 
placement in the position of molecules between t and t + At given by (1.7). One key 
assumption that is made in equation ( |2.1[ ) is that molecules that are absorbed into Sic 
from fijvf are absorbed into the closest compartment Co.j to their position calculated 
at t + At (see the limits of integration over the variable y in (2.1)). It can also be 
shown [6j [IT] that the contributions of the two ways a molecule may be absorbed by 



the interfacial compartment (see cases (a)-(b) in Section 1.3) to the total number 
of absorbed molecules in a time step are equal giving rise to the factor of 2 in the 



integral term on the right hand side of equation (2.1 ) and the absorption component 



of the integral on the right hand side of equation (2.2) 
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We may further simplify equations ( |2.1| )-( |2~2"| by using the symmetry of the do- 
main in the y-direction. Firstly, since $j (where now j refers to compartments C$j 
on the interface) is dependent only on the morphology of Cqj and its relative position 
to /, we expect <&., to be independent of j and we therefore denote 

$, = $. (2.3) 

Secondly, symmetry in the y-direction also allows us to make the conclusion that 

fj(x,y) = f {x,y - jh). (2.4) 

where the index of fj refers only to interfacial compartments Cq,j. 

In the vicinity of x = there is a boundary layer of width 0(y/A~t) [51 ITT] . 



We rescale equations (2.1|-(2.2| using the dimensionlcss boundary layer coordinate 
X = £\/DAt. We also denote the probability density and placement function in 
this boundary layer region by Pi n (£,y,t) and fmj(£,y) = VDAt fj(£y/DAt, y). The 
rescaling of fj is necessary to satisfy the normalisation condition 



/oo />oo 
/ fj(x,y)dxdy = l 
-oo JO 



since (as we will see) /,• vanishes outside of the boundary layer. Thus, using (2.3 1 and 



(2.4), in the boundary layer coordinates equations (2.1|-(2.2| become 



p 0<n (t + At) =(l- ^2 )po,n(*) + Ta(pO,n-l(*) +P0,«+l(*) + Pl,n(*)) 

9\T/ p7ih+h/2 poo /*oo poo 

+ 77 dy dy dd d£ p in (£, y, t)K(t + £)K(y-y), (2.5) 

niV Jnh~h/2 J-oo JO JO 



/OO /"OO 

dy / df pfrfa V, t) [k (e - €) + (1 - 2*) K (£ + £)] K(y-y) 
-oo JO 

+ >/5At*53/ in ,o(e > y-j/»)poj(t) 1 (2.6) 

where A = h/s/DAt and «;(£) = \/DAl K (y/DM £) = (4tt)-^ 2 exp(-£ 2 /4). Let us 
denote p(y, t) = P(0, y, t) and Pa(y, i) = P x (0, y, t) where P(x, y, t) is the distribution 
which the TRM approximates (for x € M. and j/gR). In order for the models to join 
smoothly at the interface / we require on the compartment-based side that 

Po,n(*) ~ -P(0, nh, t) = p(nh, t), 

Pi, n (t) ~ P(-h, nh, t) = p(n/i, t) - hp x (nh, t) + 0(h 2 ), 
Po, n +i(t)~P(0,nh + h,t) = p(nh + h,t), 
Po,n-i(t) ~ P(Q,nh- h,t) =p(nh- h,t), 

while, for the molecular-based side, we require no rapid variation in the boundary 
layer, so that 

p in (£, y, t) ~ p( y , i) + v / DAt(C + C x )fe(y, t) + . . . , 
p in (C, y, i + At) ~ p(y, t + At) + v A DAt(C + C x )p x {y, t) + . . . , 
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where we have allowed for a small shift C x in which the molecular-based region "sees" 
the interface [TT]. Similarly 



K{y-y) p(y, t) Ay = p(y, t) + O(At). 



Substituting the expansions (2.7)-(2.8) into equations (2.51 and (2.6 1, using h 
^/DAt as h — > and At — » 0, we obtain, 



2* ] 



= -j^p(nh,t) + -^hp x (nh,t) 
2* 



1 
AT 






= -*erfcl |lp(v,t) 

+ T^At f^— |^e-« 2 / 4 - C^erfc f|) + (* - l)£erfc (| 
+ VDAt^J2 finAtv - jh)p(jh,t) + 0(h 2 ). 



Equating coefficients of p(nh,i) and p x (nh,t) in (2.9| gives 



$ = 



2*A 



1 = * 1 



^O t 



From the coefficient of p x in (2.10) we see that 



* = 1, and C T = 0. 



(2.9) 

Px(y,t) 

(2.10) 

(2.11) 
(2.12) 



Consequently, equation (2.11) implies that $ satisfies ( |1.11[ ). Using (2.12) and ( 2.10[ ) 
we find that / in has to satisfy 

^li)p(y,t) = ^y2f in ,o(^y-jh)p(jh,t) + 0(h 2 ). 



Writing 



this becomes 



/in ) o(e ! y) = ^erfc(|)F(,i. 



P(y,t) = hJ2F(y-Jh)p(jh,t) + 0(h 2 ). 



(2.13) 



This is a standard interpolation problem. The simplest weight function which gives 
0(h 2 ) accuracy is the triangle function 



,-,,„ = < K 1 "^ 1 



(2-14) 



0. 



otherwise. 
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Thus 




Mx,y) = { hyADM^ivmmJV-i)' ^- h <y <h > (2 .i 5 ) 

otherwise. 



It is important to note that our expansions (2.7) assume that the probability pij of 
being in compartment C,j should be a continuous extension of p(x,y) evaluated at 
(—ih, jh) which is in the center of the right side of the compartments. If the expansions 



(2.7) were taken in the center of the compartments (evaluated at (—ih — h/2,jh)) 
the result would be that C\ = —A/2 and thus a shift in the continuous expected 
probability density curve over the interface is seen resulting in an apparent error of 
hp x /2 + 0(h 2 ) as At — > on the interface [TT]. This error can therefore be reduced 
by refining the compartments near the interface (remembering that this might mean 
also reducing At such that DAt ~ h 2 ). 

The derivation above was conducted under the assumption DAt ~ h 2 . Let us now 



assume DAt < h 2 instead. Then A = hj\jDAt and $ (given by (1.11) as 2A/ ' ypK) 
are no longer of order 1, and the above derivation may fail. In particular, we retain 
terms of order h while ignoring terms of order $ft 2 . If we are interested in very small 
At the errors associated with this mismatch between At and h may be reduced by 
replacing F(y) by the step function 

, , for - h/2 < y < h/2, 

'run = { /< ' ' (2.16) 

otherwise, 



which gives O(h) accuracy to (2.13). This function is also easier to implement in the 



case of corners which will be discussed in the next section. Using (2.16) instead of 



(2.14), equation (2.15) reads as follows 




(2.17) 



These issues will be discussed further in Section [3l 



2.2. Interface corners. The analysis in the previous section does not extend 
trivially to the case when there are corners in the interface I. Indeed, when the 
interface is not perfectly flat the previous section is invalid. However, since we are 
restricting our analysis to regular cubic lattices in flc, the interface must be made 
up of a series of straight edges connected at right angles. Let us consider the com- 
partments Cij assigned to the region ((i — l)h, ih) x ((j — l)h,jh), where the indices 
(i,j) € Z x Z \ N x N (i.e. they fill the complement of the positive quadrant). The 
considered geometry is represented graphically in Figure |2.2[a). 



Molecules cannot leave compartment Co,i using (2.14). This is because there is 
a non-zero probability that the molecule crosses diagonally to the region of C^o and 
vice versa. This diagonal movement is prohibited by the rules of the compartment 
region but is allowed in the molecular region. Furthermore, molecules in Qm may 
move into Co,o during one time step. Typically, molecules in compartments will not 
be able to move diagonally out of Co,o- Whilst it may be possible to make diagonal 
motion an exception for Cq.o, atypical functions placing molecules into VLm from Co,cb 
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Fig. 2.2. (a) Illustration of corner geometry in two-dimensions. Compartment indices are 
denoted as in Section |2.2| (b) The corner geometry used in illustrative numerical simulations in 
Section I3.2I 



Co,i an d Ci.o must be determined and also rules for how molecules in £Im close to the 
corner migrate diagonally, left or down into £lc must be determined. It is important 
to note that the added complexity to these corner compartments is not trivial like the 
case when the compartments form a straight interface. If we attempted to generalize 



(2.14), then complex distributions would have to be sampled from to place molecules 
from Co.i and C\$ into Qm and several tests would have to be performed on molecules 
in Q,m (to see if they meet criteria for diagonal migration into flc, or to determine 
if they move down or left over the interface). Given this complication, we find it 



reasonable to use (2.16) instead of (2.14) which means that we use (2.17) instead of 
(2.15). Distribution (2.17) does not allow for molecules to leak between Cq,i and Ci : o 



but at the cost of accuracy in the local region around the corner. 

The treatment of molecules at the corner is therefore described by the following 



rules. Molecules in compartments Cq.i and C\$ have the propensities given by (1.10) 
namely 



D , 



*pM,o, 



where A/o,i (resp. A/i,o) is the number of molecules in the compartment Cq.i (resp. 
C\fi). Molecules from Co,o ma y not enter J7m directly. Molecules from Q,m that land 
in Co,o or satisfy condition (b) in Section |1.3.2 for both parts {x = 0} and {y = 0} of 
the interface / are placed at random (with probability 1/2) in Cx,o ° r £-0,1- Molecules 
in the compartment Cq.i migrating to Qm are placed according to the distribution 
(2.17). Molecules in the compartment C\$ migrating to Qm are placed according to 



the distribution 



fo(x,y) 




erfc 



(_V 

VV4DA* 



for - h/2 < x < h/2, 
otherwise, 



(2.18) 



which can be obtained from the distribution (2.17) by exchanging the variables x and 
y. Thus, in both cases, we use (1.12) perpendicular from their respective interfaces 



and the step distribution (2.16) tangentially along each respective interface. 

Since DAt < h 2 molecules that are in compartments Ci.o and Cq.j (i, j > 2) are 
not significantly affected by the corner within one time step. We therefore use the 
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probability distribution (2.171 for these compartments. Thus we use the distribution 



(2.16) for molecules leaving the corner compartments, and the triangle distribution 



(2.14) for molecules leaving all other boundary compartments. 



Finally, we note that we also tested the alternative of using a one-sided triangle 
distribution for the corner compartments C\_q or Co.i- It produced results which were 



indistinguishable numerically from the distribution (2.16), although this one-sided 



distribution includes a small bias of particles away from the corner. 

2.3. Parameters for a flat interface in N dimensions. The derivation pre- 
sented in Section [2~T] can be used for flat interfaces in arbitrary dimensions. One can 
show that the parameters Qj = $ and W are independent of the number of dimensions 
N of the domain for a flat interface with a regular cubic compartment arrangement 
in flc (and therefore take the values derived in Section 2.1). Furthermore one can 



show that the distribution /o(x) for placing molecules in VLm is the product of N 
distributions separating the coordinates 



A 



fo(xi,x 2 ,...,x N ) = F J -(xi)l[FU 



(2.19) 



i=1 



Here, the distribution F ± (xi) is for the coordinate X\ perpendicular to the interface 
/ = {xi = 0}. It is given (see equations (1.12) and (2.15)) by 



F ± (x 1 ) = 



-erfc 






j'i 



4DAt \V4DAi 



The remaining N — 1 identical distributions for each coordin ate x 



(2.20) 



,N, 



tangential to the interface I, are denoted as F"(xi) in equation (2.19). If the origin is 
placed at the center of c ompartment in question, then one c an fol low the deri vation 
presented in Section 2.1 to obtain F"(xi) — F(xi) given by (2.14). Equation (2.19) 



indicates that tangential coordinates should be chosen independently from each other. 

3. Illustrative numerical examples. In this section we present simple diffu- 
sion simulations using the TRM to demonstrate its accuracy and convergence. We 
will compare results computed by distributions (2.15) and (2.17). 



3.1. Straight interface. We will consider a diffusing molecule which starts at 
the origin at (dimensionless) time t = and diffuse with (dimensionless) diffusion 
constant D = 1 in the semi-infinite two-dimensional domain tt — (0, oo) x (0, oo). 
Boundary dft will be considered reflective. 

The probability distribution, P(x,y,t), to find the molecule at time t given its 
initial position at the origin evolves according to the partial differential equation 
(PDE) 



dP 
~dt 



d 2 P d 2 P 



dx 2 dy 



2 ' 



(3.1) 



with the initial condition P(x, y, 0) = 5{x, y) where 8 is the Dirac delta function. 
Using no-flux (reflective) boundary condition (VP • n 
solved as 



0), equation (3.1) can be 



P(x,y,t) = — exp 

wt 



-(x 2 + y 2 ) 
At 



(3.2) 
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Fig. 3.1. Probability distribution at time t = 0.5 estimated using No = 2 X 10 J realisations 
of the TRM method for the domain partition pOt . (a) The expected distribution found by ( |3- 2[ ) . 
(b) TRM simulation with compartment size h = 0.05 and At = 0.0004. (c) TRM simulation with 
compartment size h = 0.1 and At = 0.0016. (d) TRM simulation with compartment size h = 0.25 
and At = 0.01. Clc can be seen in (b)-(d) o n the right and Sin to the left of the white solid line. 
These simulations were done using sampling 12.151. 



We shall simulate this diffusion process stochastically using the TRM where Q, is 
divided into 



Q M = (0,0.5) x (0,oo) 



and 



Q c = (0.5, oo) x (0,oo). 



(3.3) 



We run Nq = 2 X 10 5 realizations of the TRM for the molecule starting at the origin 
and note its state (either in the molecular regime or compartment regime) at each time 
step until t — 1. We simulate the compartment regime using a square lattice with non- 
dimensional compartment spacings h = 0.05, h = 0.1 and h = 0.25. Since our analysis 
uses the assumption DAt ~ h 2 , we use the time steps At = 0.0004, At = 0.0016 and 
At = 0.01 respectively with D = 1 such that h/y/DAt = 2.5 - 0(1). 

At t = 0.5 the molecule positions (for each realisation) are binned according to 
their compartment (or in the case of the molecular regime, counted in bins of area 
h 2 ) and a plot of these bin copy numbers divided by N^h 2 is produced to show the 
approximate probability density that is generated by the TRM. These probability 



densities arc shown for comparison against the exact solution (3.2) in Figure 3.1 for 



each compartment size h. The distributions generated using the TRM match well 
with the expected distribution and appear to be more accurate as h is decreased. 
To better visualize the accuracy of the TRM, we define the error function 



Error(t) 



G 



TRM 



(t) 



m 



P(x,y,t) dx dy, 



(3.4) 



where Ctrm (t) is the number of realisations of the TRM which have the molecule po- 
sitioned in $7<7 at time t. Therefore, the fraction Ctrm^/Nq is the approximation of 
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Fig. 3.2. (a) The error of the TRM defined by pTI} for h = 0.05 and At = 0.0004 (blue line), 
h = 0.1 and At = 0.0016 (red line) and h = 0.25 and At = 0.01 (black line). TRM simulations are 
taken with Qq = (0.5, oo) X (0, oo) and Qm = (0,0.5) X (0, oo). These simulations were done using 
sampling ||2.15|. (b) Results from panel (a) scaled by h. 



I In P( x > Vi ^) d% dy and the error ( 3.4 1 measures the accuracy of this approximation. 



In Figure |3T2[ a) we present the error (3.4) as a function of time for the h = 0.05, 
h = 0.1 and h = 0.25 simulations. There is a maximum in this error around t ~ 0.1 for 



all simulations. The predicted error on the boundary (see Section 2.1 ) is proportional 
to the net flux over the interface (hp x /2 + 0(h 2 )). We note that this flux reaches 
a maximum (according to the exact solution (3.2)) at around t — 1/24. The reason 



why the maximum of our measured error does not match up with this time is because 
the error at the interface is described by a small discontinuity in the distribution on 
the boundary, this discontinuity then diffuses into each of the subdomains. After 
t = 1/24 the discontinuity is reduced, thereby reducing this bias effect but there is 
still some time before the distributions are corrected by diffusion. It is for this reason 
that coupling at the interface correctly is so crucial. A small bias in the flow over the 
interface in one direction can lead to an avalanching effect on the distribution. Whilst 
the TRM cannot eliminate the error that is associated with changing of regime entirely, 
it does optimize the error. We also expect the error, however, to be proportional to 
h. We can see that this is the case by plotting the error ( |3.4[ ) divided by h for each 
simulation (Figure 13721(b)). Each of the three curves in Figure [3~2) (b) approximately 
overlay implying that the leading order term of the error is 0(h). This means that in 
the continuous limit h — > the error that is due to the TRM appoximately converges 
linearly with h, as expected. 

The purpose of the TRM is to match the concentration and perpendicular flux of 
molecules at the boundary in the most optimal way given small h and small At. In 
our analysis, we restricted the algorithm parameters to the case DAt ~ h 2 . In some 
applications, this might not be a preferred parameter regime because At determines 
the resolution of the microscopic region SI m and if DAt ~ h 2 then this resolution is no 
better than the compartment-based regime. One must be careful in this case since for 
DAt <C h 2 , the parameter $ = 2A/y / 7r (defined by (2.11)) may become significantly 



larger than 1. Previously, higher order terms of the form $/i 2 in (2.9)— (2.10) (which 



were ignored as "too small" ) now become dominant over terms that are of the order 
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of h. This manifests itself into an increase in the dispersion inside the boundary layer 
near the interface in the tangential direction. This effect can be seen in Figure [3~3"t c). 
The increase in the dispersion in the tangential direction may also be explained 
by following the TRM mechanism. Consider molecules in Cl^j close to the interface 
between S!m and Clc- If these molecules are sufficiently close to the interface, they 
are likely to be absorbed into a compartment in ilp. To compensate for this rapid ab- 
sorption of molecules, the compartments on the interface must return these molecules 
to the boundary layer at comparable rates. Since, as At — > 0, the boundary layer 
gets thinner, this increases the rate of resorption of molecules having just come from 
flc- If a molecule enters an interfacial compartment Coj near the boundary to the 
adjacent compartment Coj-i then this molecule effectively jumps a distance h/2 in 
the tangential direction (because its position in flc can be considered as the center of 
Cqj). This molecule is then rapidly interchanged between f2jvf and Cqj each time with 
a new tangential coordinate until the tangential coordinate falls out of line with Cj or 
the molecule diffuses away from the interface. The former of these two options occurs 
more often if At is small since the tangential coordinate is rapidly resampled until the 



0.25 chance of being sampled outside the compartment Cj is realized if we use (2.14). 



In these situations, it is better to use (2.16) for the tangential coordinate. This is 



because, (2.16) restricts the molecule initiation to the compartment from which it 



came. Furthermore, one can reduce this effect by reducing the size of h. Reducing h 
helps to reduce $. 

It is possible to show that this numerical artefact is improved, if h cannot be 



reduced, by using step function (2.16) instead of triangle function (2.14) to s ample 
molecule positions tangentially to the interface. However, if DAt ~ h 2 then (2.14) 
offers the best results. This is demonstrated in Figure |3.3| In Figure |3.3( the least 
amount of artificial dispersion is achieved if DAt ~ h 2 using the triangle function 



sampling (2.14) but the triangle function sampling introduces more severe artificial 
dispersion along the interface than step function sampling (2.16) if DAt <C h 2 . 



3.2. Interface with a corner. To demonstrate that the TRM produces good 
results when there is a corner in the interface, we present also the results of TRM 



simulations of the same problem as in Section 3.1 with subdomains redefined as follows 



Sl M = (0,0.5) x (0,0.5), ft c = (0,0.5) x (0.5, oo)U (0.5, oo) x (0,oo). (3.5) 



The details of the TRM implementation of corners are discussed in Section [272) The 
corner is oriented as in Figure [272[b), i.e. the TRM implementation presented in 



Section 2.2 for the corner orientation in Figure 2.2 'a) is adjusted (by a simple rotation) 



to the corner orientation in Figure |2.2[ b) . 

In Figure |3.4| we present the distributions that result using the TRM for the 



domain partition (|3.5|). These distributions are calculated in the same way as those 

0.5. 



distributions found in Figure |37l) The distributions are plotted, similarly, at t 
There is good agreement with the expected distribution (Figure |3~4| a)) especially for 
small h. 
h- 



In Figure 3.5 'a), we present the error (3.4) as a function of time for the 
h 



0.05, h = 0.1 and h = 0.25 simulations. To verify that the error due to the TRM 
is still O(h) when the interface has a corner in it, Figure [3~5^ b) shows the error for 
the TRM simulations scaled by h. 



In Figure 3.6 'a), we again see the artefact for small At in the compartment 



described by the region x € [0.5,0.75), y € [0.5,0.75); an unexpectedly large number 
of molecules. This is because diffusion is biased tangentially along both sides of the 
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A t « h 2 



At ~ h 3 




if) 



M 



Fig. 3.3. Probability distribution at time t = 0.5 estimated using Nq = 2 x 10 5 realisations of 
the TRM method for the domain partition ( |3 - 3[ ) - In all simulations h = 0.25 and D = 1. f2o con fee 
seen in (a)-(d) on the right and Qm to the left of the white solid line denoti ng int erface I . TRM 
simulations with (a) DM = 10~ 5 < h 2 and l |2.15f ; (b) DAt = 0.0 1 ~ /i 2 and |2.15| ; 



(c) DAi = 10" 5 < /i 2 and ( |2.17| ,- (d) DAt : 



: 0.01 ■ 



: 0.0 1 
h 2 and (|2.17k. 



interface (due to the large h and small At as we discussed in Section 3.1) and these 
molecules gather at the corner. 



4. Discussion. In this paper we presented an analysis of the TRM on regular 
lattices in dimensions larger than one. We have derived TRM parameters to simulate 
diffusing molecules migrating over interface / that separates domain flc in which 
molecules are constrained to a lattice and domain VLm whereby molecules can diffuse 
in continuous space. The TRM algorithm is presented in Table |1.2| We considered 
two cases DAt <C h 2 and DAt ~ h 2 . We showed that the parameter 3>jj defined 



in (1.10) can be chosen as in the one-dimensional case (1.11). The distribution /(x) 
for placing molecules in Qm was presented in equation (2.19). In particular, each 



tangential direction to the interface may be treated independently and in the same 



way. We presented two different approaches (2.14) and (2.16) for sampling tangential 
directions. 

In the case DAt <C h 2 , the step function approximation (2.16) is the most ap- 
propriate choice. In particular, the distribution for placing molecules in Qm is given 



by (2.17) for two-dimensional problems. In the case DAt ~ h , the triangle function 



approximation (2.14) is the most appropriate choice, i.e. the distribution for placing 
molecules in J7j\/ is given by (2.15 1. Moreover, Figure 3.3 demonstrates that the overa ll 
best results (in terms of accuracy) can be obtained in the case DAt ~ h 2 and (2.14). 
However, imposing the condition DAt ~ h 2 over the whole domain ile might lead to 
computationally intensive simulations. A natural solution to this problem would be 
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Fig. 3.4. Probability distribution at time t = 0.5 estimated using TVq = 2 X 10 5 realisations of 
the TRM method for the domain partition ( |3- 5[ ) - (a) The expected distribution found by the exact 
solution |3.2[ l. (b) TRM simulation with compartment size h = 0.05 and At = 0.0004. (c) TRM 
simulation with compartment size h = 0.1 and At = 0.0016. (d) TRM simulation with compartment 
size h = 0.25 and At = 0.01. Vic can be seen in (b)-(d) outside of VIm which is boxed by th e white 
solid line (interface I) about the origin. These simulations were done using sampling l |2,17| l- J2,18[ l 
near the corners according to the approach described in Section\2.2\ 



to use unstructured meshes j5j where compartments can be of different sizes. Then 
we could impose the condition DAt ~ h 2 for compartments close the interface I (to 
maximise accuracy) and the condition DAt <C h 2 for compartments further from the 
interface (to maximise efficiency). 

The last decade has seen a number of different algorithms appear in the scientific 
literature with the purpose of coupling two domains in space with different modelling 
techniques for reaction-diffusion processes similar to that of the TRM. Some of these 
algorithms aim to couple mesoscopic stochastic simulations on a lattice (compartment- 
based model) with a deterministic PDE-based mean-field description [TJ [131 1301 HE] • 
Most of these algorithms include the use of an overlap region and use this overlap 
region to calculate the flux of molecules (and other conserved physical quantities) 
that are flowing between the two regimes. This flux is calculated by particle counting 
and matched to the flux condition which forms the boundary condition of the deter- 
ministic macroscopic approach [TJ [T31 [301 I2S] • New particles are created either in a 
new compartment corresponding to a lattice-point used in the PDE solver that is im- 
plemented on the deterministic side of the interface ( "handshaking region" ) p] or in 
the closest compartment to the position at which the deterministic region terminates 
[131 130] . Whilst these approaches often have a region of overlap where the description 
of particles is transitioning from deterministic to compartment-based, these regions 
are typically thin and transition between the descriptions is clearly defined in space. 
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Fig. 3.5. (a) The error of the TRM defined by §SM for h = 0.05 and At = 0.0004 (blue line), 
h = 0.1 and At = 0.0016 (red line) and h = 0.25 and At = 0.01 (black line). TRM simulations are 
taken with VLm and Vic defined in ||3.5|. (b) Results from (a) scaled by h. 
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Fig. 3.6. Probability distribution at time t = 0.5 estimated using Nq = 2 x 10 5 realisations of 
the TRM method for the domain partition ( |3.5| . In both simulations h = 0.25 and D = 1. Vic can 
be seen outside ofVlM which is boxed by the white solid line denoting interface I. TRM simulations 
with (a) DAt = 10 — 5 <C h 2 ; (b) DAt = 0.01 ~ h? . These simulations were done using sampling 
||2.17|-J2,18| near the corners according to the approach described in Section\2.2\ 



Ferm et al. [10] presented a new technique for connecting deterministic regions with 
mesoscopic regions using a smoother coupling technique. Their technique involves in- 
terfacing a region of deterministic description with a stochastic description modelled 
using a tau leaping algorithm. The tau leaping algorithm is a quick compartment- 
based algorithm that does not account for all events as they occur but rather updates 
the system at discrete moments in time r [18] . This stochastic simulation technique 
is fast but is only accurate if the concentration of molecules is large enough that in- 
dividual events on smaller time scales produce small perturbations that do not affect 
the simulation significantly. This tau leaping algorithm is then interfaced with a more 
accurate compartment-based algorithm as the concentration drops. 
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Geyer et al. |15[ I19| developed a method for coupling a deterministic numeri- 
cal solution to a reaction-diffusion PDE with a Brownian dynamics molecular-based 
algorithm. In a similar way to the previous algorithms which couple deterministic 
PDEs and compartment-based algorithms the flux is determined at the interface. 
Molecules are initiated in the molecular-based algorithm in one-dimension using the 



same distribution that is stated in this manuscript (1.12 1 (see (11) in [E]). Indeed 
this boundary condition is generated by considering a "compartment" with a partic- 
ular expected number of molecules. Absorption of molecules at the interface is done 



without consideration of condition (b) in Section 1.3.2 Since this condition is neces- 
sary in the TRM to generate the correct expected flux on both sides, the matching of 
the flux is instead forced manually by the algorithm. Unlike the algorithm presented 
in [TS] the TRM couples two stochastic simulations. In the case of the TRM, the ex- 
pected flux cannot be calculated without averaging over time. We present the TRM 
specifically with rules that govern each molecule as they cross the interface rather 
than try to introduce conditions that depend on the simulation itself. This is a more 
natural methodology since it should not be the case that an individual molecule is 
influenced by the net flux but rather be treated independently from each other. 

Franz et al. [2] recently developed a technique for coupling deterministic PDEs 
with microscopic molecular-based simulations. Whilst this technique is different from 
the TRM as it interfaces a deterministic PDE-based description with a stochastic 
one, it deserves special mention because unlike other techniques that couple deter- 
ministic systems with stochastic simulations it does not use an averaging technique to 
determine the flux over the interface. Rather, it treats each particle as an individual 
as it crosses from the molecular-based simulation and is added to the deterministic 
description in a probabilistic sense, where the probability for finding the particle is 
known from the simulation. Transversely, molecules may migrate back when they 
are sampled from the probability distribution which is proportional to the continuous 
distribution of particles. 

An algorithm similar to the TRM was recently published by Klann et al. [22] . In 
this algorithm, compartment-based and molecular-based regions are coupled. Klann 
et al. [22] present no comparison with the TRM and it is not clear how exactly they 
couple different regimes. However, several differences between their approach and the 
TRM can be identified. 

The analysis presented in [TT] and in Section [2] reveals that some steps of the 
Klann et al. hybrid method [22] are not optimal. For example, in [22], a molecule 
migrates from the compartment-based algorithm with a propensity which is natural 
for the lattice. This means that they have taken a value of <& = 1 instead of $ given 



by (|1.11[). Instead of placing the molecule with a distribution given by (2.151 into the 



the molecular-based domain they place the molecule within the compartment that 
corresponds with a natural extension of the compartment-based approach. Further- 
more, molecules are transferred back into compartment-based molecules without the 



condition (b) in Section 1.3.2 (this corresponds to a value of \P = 1/2 instead of \1/ 
given by (1.11 1). It is our experience that somewhat ad hoc or heuristic coupling 
of this type between compartment-based and molecular-based regions, despite being 
simpler to implement, may lead to heavy biasing of molecules especially in the case 
when the expected net flux over the interface is high. This biasing also depends on 
the relationship between At and h. 
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